library(reshape2)
library(ggplot2)
theme_set(theme_classic() +
		theme(text=element_text('serif')) +
		theme(strip.background=element_blank()) +
		theme(strip.text.x=element_text(face='bold',size=12)))

drawmodel <- function(ode_fn, ssa_nbr, ssa_fn, export_fn, labels)
{
	X.ode <- ode_fn()
	X <- cbind(melt(as.data.frame(X.ode), id='time'), type='ode', group=0, insample=FALSE)

	library(parallel)
	cl <- makeCluster(detectCores()/2)
	clusterExport(cl=cl, varlist=c("ssa_fn","model.solve","model.euler","model.tau_leaping","model.maruyama",export_fn), envir=environment())
	nsim <- 400  #80  #20
	X.ssa <- parLapply(cl, 1:nsim, function(i) ssa_fn())
	stopCluster(cl)

	j <- sample(100)
	
	for(i in 1:length(X.ssa)) {
		if(X.ssa[[i]][100/0.01,'V'] != 0)
			X <- rbind(X, cbind(melt(as.data.frame(X.ssa[[i]]), id='time'), type='ssa.gt0', group=i, insample=any(i==j)))
		else
			X <- rbind(X, cbind(melt(as.data.frame(X.ssa[[i]]), id='time'), type='ssa', group=i, insample=any(i==j)))
	}

	#X.avg$value <- loess(X.avg$value ~ X.avg$time+X.avg$variable)$fitted
	#X.avg$value <- loess(value ~ variable, X.avg, subset=(variable=='T'))$fitted
	
	X.avg <- aggregate(value ~ time+variable, subset(X,type=='ssa'), mean)
	X <- rbind(X, cbind(X.avg, type='avg', group=length(X.ssa)+1, insample=FALSE))

	X.avg <- aggregate(value ~ time+variable, subset(X,type=='ssa.gt0'), mean)
	X <- rbind(X, cbind(X.avg, type='avg.gt0', group=length(X.ssa)+2, insample=FALSE))

	X.avg <- aggregate(value ~ time+variable, subset(X,type=='ssa'&insample==TRUE), mean)
	X <- rbind(X, cbind(X.avg, type='avg1', group=length(X.ssa)+3, insample=FALSE))
	X.avg <- aggregate(value ~ time+variable, subset(X,type=='ssa.gt0'&insample==TRUE), mean)
	X <- rbind(X, cbind(X.avg, type='avg1.gt0', group=length(X.ssa)+4, insample=FALSE))
	
	
	
	X[X$type=='ssa.gt0',]$type <- 'ssa'

	X$variable <- factor(X$variable,
		levels=levels(X$variable),
		labels=labels)

	library(plyr)  # for subset, to limit yaxis

	print(ggplot(X, aes(time,value,group=group,alpha=type,size=type,color=type)) +
			#geom_line(subset=.(variable=='Target cells'|(value<5&variable=='Infected cells')|(value<15000&variable=='Free virus'))) +
			geom_line() +
			facet_wrap(~ variable, scales='free_y') +
			expand_limits(y=0) +
			scale_alpha_manual(labels=c('ode','ssa','avg','avg.gt0','avg1','avg1.gt0'), values=c(1,min(1,2.5/nsim),1,1,1,1)) +
			scale_size_manual(labels=c('ode','ssa','avg','avg.gt0','avg1','avg1.gt0'), values=c(0.5,0.3,0.3,0.3,0.3,0.3)) +
			scale_color_manual(labels=c('ode','ssa','avg','avg.gt0','avg1','avg1.gt0'), values=c('black','black','red','blue','darkred','darkblue')) +
			xlab('Time units') +
			ylab('Cell/viral units') +
			theme(legend.position='none'))
}
